---
title: "Vignette"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

require(memisc)
require(stringi)
library(stargazer)
library(coefplot)
library(interplot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(stargazer)
library(stringr)
library(tidyverse)
library(ggplot2)
library(psych)
library(margins)
library(dplyr)
library(modelsummary)
library(makedummies)
library(ggeffects)
library(estimatr)
library(ggpubr)
library(BalanceR)

library(haven)
library(summarytools)
library(tidyverse)
library(ggplot2)

dat <- read_csv("data/Iraq2024.csv")
dat <- dat[-1,]

dat$Q17a[is.na(dat$Q17a)] <- 0
dat$Q17b[is.na(dat$Q17b)] <- 0
dat$Q17c[is.na(dat$Q17c)] <- 0

dat <- dat |>
  mutate(Support = Q17a + Q17b + Q17c)
dat$Support[dat$Support == 0] <- NA

dat <- dat |>
  mutate(Treatment = case_when(Q17a >= 1 ~ "Severer",
                               Q17b >= 1 ~ "Scrapped",
                               Q17c >= 1 ~ "Tolerance"))
dat$Treatment <- as.factor(dat$Treatment) 

dat <- dat |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5)

dat <- dat |>
  mutate(sunni = ifelse(D6 == 1, 1, 0),
         shia = ifelse(D6 == 2, 1, 0),
         christian = ifelse(D6 == 3, 1, 0),
         kurd = ifelse(D6 == 4, 1, 0),
         others = ifelse(D6 == 5, 1, 0))

dat <- dat |>
  mutate(Sect = dplyr::recode(D6, 
                       "1" = "Sunni",
                       "2" = "Shia",
                       "3" = "Christian",
                       "4" = "Kurd",
                       "5" = "others"))

dat <- dat |>
  mutate(Islamism = Q3_4)
dat$Islamism[dat$Islamism == 6] <- NA
dat$Islamism <- 6 - dat$Islamism

dat <- dat |> 
  mutate(islamism = case_when(Q3_4 <= 2 ~ "High",
                              Q3_4 == 4 | Q3_4 == 5 | Q3_4 == 3 | Q3_4 == 6 ~ "Low"))

dat <- dat |>
  mutate(Liberalism = Q3_5)
dat$Liberalism[dat$Liberalism == 6] <- NA
dat$Liberalism <- 6 - dat$Liberalism

# Liberal_high_education
dat <- dat |> 
  mutate(education = case_when(D3 <= 2 ~ "Low",
                               D3 == 3 | D3 == 4 ~ "Middle",
                               D3 >= 5 ~ "High"))
dat <- dat |>
  mutate(liberalism = case_when(Q3_5 <= 2 ~ "High",
                                Q3_5 == 4 | Q3_5 == 5 ~ "Low",
                                Q3_5 == 3 | Q3_5 == 6 ~ "other"))

dat <- dat |>
  mutate(Liberal_high_education = ifelse(liberalism == "High" & education == "High", "High-educated_liberal", "Other"))


dat <- dat |> 
  mutate(Provincial_capital = case_when(D8 == 1 ~ "capital",
                                        D8 == 2 ~ "rural"))
dat <- dat |> 
  mutate(Shia_party_supporter = (Q4_1 + Q4_2 + Q4_3)/3) |>
  mutate(Shia_party_supporter = ifelse(Shia_party_supporter >= 6, 1, 0))
```

# balance check
```{r}
balance <- BalanceR(data = dat, 
                    group = Treatment,
                    cov = c(Gender, Age, Education, Income))
print(balance, digits = 3)
summary(balance, digits = 5)
plot(balance, point.size = 5, text.size = 18)
```

# boxplot
```{r}
dat |> 
  ggplot() +
  geom_boxplot(aes(x = Treatment, y = Support)) +
  scale_x_discrete(na.translate = FALSE)
```

# regression
```{r}
reg01 <- lm(formula = Support ~ Treatment,
            data = dat)
summary(reg01)

reg02 <- lm(formula = Support ~ Treatment +
              Gender + Age + Education + Income,
            data = dat)
summary(reg02)

reg03 <- lm(formula = Support ~ Treatment +
              Sect +
              Gender + Age + Education + Income,
            data = dat)
summary(reg03)

reg04 <- lm(formula = Support ~ Treatment +
              Islamism +
              Gender + Age + Education + Income,
            data = dat)
summary(reg04)

reg05 <- lm(formula = Support ~ Treatment +
              Provincial_capital +
              Gender + Age + Education + Income,
            data = dat)
summary(reg05)

reg06 <- lm(formula = Support ~ Treatment +
              Shia_party_supporter +
              Gender + Age + Education + Income,
            data = dat)
summary(reg06)
```

# regression: interaction terms
```{r}
reg07 <- lm(formula = Support ~ Treatment * Sect +
              Gender + Age + Education + Income,
            data = dat)
summary(reg07)

reg08 <- lm(formula = Support ~ Treatment * Islamism +
              Gender + Age + Education + Income,
            data = dat)
summary(reg08)

reg09 <- lm(formula = Support ~ Treatment * Provincial_capital +
              Gender + Age + Education + Income,
            data = dat)
summary(reg09)

reg10 <- lm(formula = Support ~ Treatment * Shia_party_supporter +
              Gender + Age + Education + Income,
            data = dat)
summary(reg10)
```

# regression: chategorical variables of Islamism and Liberalism
```{r}
reg11 <- lm(formula = Support ~ Treatment +
              relevel(factor(islamism), ref = "Low") +
              Gender + Age + Education + Income,
            data = dat)
summary(reg11)

reg12 <- lm(formula = Support ~ Treatment * relevel(factor(islamism), ref = "Low") +
              Gender + Age + Education + Income,
            data = dat)
summary(reg12)

reg13 <- lm(formula = Support ~ Treatment + 
              Liberalism +
              Gender + Age + Education + Income,
            data = dat)
summary(reg13)

reg14 <- lm(formula = Support ~ Treatment * Liberalism +
              Gender + Age + Education + Income,
            data = dat)
summary(reg14)

reg15 <- lm(formula = Support ~ Treatment * Liberal_high_education +
              Gender + Age + Education + Income,
            data = dat)
summary(reg15)
```

# regression table
```{r}
reg_table <- mtable("Model 1" = reg01, "Model 2" = reg02, 
                    "Model 3" = reg08, "Model 4" = reg14,
               summary.stats=c("adj. R-squared","F","p", "N"))
print(reg_table)
write_html(reg_table, file = "result/result_vignett.html")

robust_table <- mtable("Model 3" = reg12, "Model 4" = reg15, 
               summary.stats=c("adj. R-squared","F","p", "N"))
print(robust_table)
write_html(robust_table, file = "result/result_vignett_robust.html")
```

# marginal effect
```{r}
plot_model(reg01, type = "pred", terms = c("Treatment"),
           axis.labels = "Treatment", title = "")

plot_model(reg02, type = "pred", terms = c("Treatment", "Gender"),
           axis.labels = "Treatment", title = "")
plot_model(reg02, type = "pred", terms = c("Treatment", "Education"),
           axis.labels = "Treatment", title = "")
plot_model(reg02, type = "pred", terms = c("Treatment", "Income"),
           axis.labels = "Treatment", title = "")
plot_model(reg02, type = "pred", terms = c("Treatment", "Gender", "Education"),
           axis.labels = "Treatment", title = "")
plot_model(reg02, type = "pred", terms = c("Treatment", "Gender", "Income"),
           axis.labels = "Treatment", title = "")

plot_model(reg08, type = "pred", terms = c("Treatment", "Islamism"),
           axis.labels = "Treatment", title = "")
plot_model(reg08, type = "pred", terms = c("Islamism", "Treatment"),
           axis.labels = "Treatment", title = "")

plot_model(reg12, type = "pred", terms = c("Treatment", "islamism"),
           axis.labels = "Treatment", title = "")

plot_model(reg14, type = "pred", terms = c("Treatment", "Liberalism"),
           axis.labels = "Treatment", title = "")
```

# prots
```{r}
p1 <- plot_model(reg01, type = "pred", terms = c("Treatment"),
           axis.labels = "Treatment", title = "")

p2 <- plot_model(reg08, type = "pred", terms = c("Treatment", "Islamism"),
           axis.labels = "Treatment", title = "")

p3 <- plot_model(reg08, type = "pred", terms = c("Islamism", "Treatment"),
           axis.labels = "Treatment", title = "")

p4 <- plot_model(reg12, type = "pred", terms = c("Treatment", "islamism"),
           axis.labels = "Treatment", title = "")

ggarrange(p1, 
          p2 + theme(axis.title.y = element_blank()), 
          nrow = 1, ncol = 2, align = "hv")

ggarrange(p1, 
          p3 + theme(axis.title.y = element_blank()), 
          nrow = 1, ncol = 2, align = "hv")

ggarrange(p1, 
          p4 + theme(axis.title.y = element_blank()), 
          nrow = 1, ncol = 2, align = "hv")
```

